/** ANAYLSIS_05.C
 *  Analyse Macros zum F-Praktikumsversuch 05: Neutronennachweis
 * 
 *  date:   25-Nov-11
 *  author: Ch.Weinsheimer 
 */

#include <iostream>
#include <stdio.h>

/**
 * User Defined Functions
 */

/// Define gauss function for fitting
double gaussian(double* vars, double* pars) {
    return (pars[3]+pars[0]*TMath::Gaus(vars[0],pars[1],pars[2]));
}

/// Define linear function for fitting
double lin(double* vars, double* pars) {
    return (pars[0]*vars[0]+pars[1]);
}

/// Read data from file and plot a histogram
void histo(char* filename)
{
    int const N = 4096;  

    // File Opening Reading and Closing
    FILE* file = fopen(filename, "r");
    char dummy[12];
    // Lese Daten
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);

    // Create 1-D Histogram
    TH1F* Hist = gROOT->FindObject("Hist");
    if(!Hist) Hist = new TH1F("Hist", "Histogramm", N, 0, N);
    else Hist->Reset("M");

    // Fill histogramm
    for(Int_t i=0; i < N; i++){
      Hist->Fill(channel[i],value[i]);
    }

    Hist->Rebin(2);

    // Setze Styles
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(0);

    // Setze Beschriftungen und Styles
    Hist->SetTitle("ADC Spectrum");
    Hist->GetXaxis()->SetTitle("channel");
    Hist->GetYaxis()->SetTitle("counts");

    // Erzeuge das Fenster
    TCanvas* Canva = gROOT->Find Object("Canva");
    if(!Canva) Canva = new TCanvas("Canva", "Versuch_05", 1024, 768);
    else Canva->Clear();
    Canva->cd(1); // Waehle das Feld an

    // Male Graph
    Hist->Draw("PL");
}



/// Read data from file and plot a histogram
void fitPeaks(char* filename,int const npeaks)
{
    int const N = 4096;  

    // File Opening Reading and Closing
    FILE* file = fopen(filename, "r");
    char dummy[12];
    // Lese Daten
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);

    // Create 1-D Histogram
    TH1F* fHist = gROOT->FindObject("fHist");
    if(!fHist) fHist = new TH1F("fHist", "Fitted Histogramm", N, 0, N);
    else fHist->Reset("M");

    // Fill histogramm
    for(Int_t i=0; i < N; i++){
      fHist->Fill(channel[i],value[i]);
    }

    fHist->Rebin(2);


    // Declare Fit-Parameters and Fit-Funtions
    double mean[npeaks], stdv[npeaks], norm[npeaks], offset[npeaks], min[npeaks], max[npeaks];
    
    TF1 fit[npeaks];

    for (unsigned int j = 0; j < npeaks; j++) {
        cout << "Setting Parameters for Peak " << j + 1 << ": " << endl << "******************************" << endl;
        cout << "Position   : "; cin >> mean[j];
        cout << "Lower Bound: "; cin >> min[j];
        cout << "Upper Bound: "; cin >> max[j];
        
        char name[2];
        sprintf( name,"%d", j);
        fit[j] = new TF1(name,gaussian,min[j],max[j],4);
        fit[j].SetParNames("Norm", "Mean", "StdDev", "Offset");
        fit[j].SetParameters(100,mean[j],100,1);
        fit[j].SetLineColor(kRed);

        fHist->Fit(&fit[j],"R0");
//        TVirtualFitter *fitter = TVirtualFitter::GetFitter();
//        fitter->PrintResults(2,0.);
        cout << "Chi-Square: " << fit[j].GetChisquare() << "/" << fit[j].GetNDF() << endl;
        
    }


    // Erzeuge das Fenster
    gROOT->SetStyle("Plain");
    TCanvas* fCanva = gROOT->Find Object("fCanva");
    if(!fCanva) fCanva = new TCanvas("fCanva", "Fitted Versuch_05", 1024, 768);
    else fCanva->Clear();
    fCanva->cd(1); // Waehle das Feld an

    // Setze Styles
    gStyle->SetOptStat(0);


    // Setze Beschriftungen und Styles
    fHist->SetTitle("Fitted ADC Spectrum");
    fHist->GetXaxis()->SetTitle("channel");
    fHist->GetYaxis()->SetTitle("counts");
    


    // Male Graph
    fHist->Draw("PL");
    for (j=0; j<npeaks; j++) {
        fit[j].DrawClone("Same");
    }
}


/// Read data from file and perform a time calibration
void timecal(char* filename)
{
    int const N = 4096;  
    Float_t peaks[12] = {702.0, 864.0, 1073.5, 1279.0, 1491.0, 1666.0, 1883.0, 2090.5, 2304.5, 2482.5, 2697.0, 2903.5};
    Float_t dpeaks[12] = {1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5};
    Float_t times[12] = {14.0, 18.0, 22.0, 26.0, 30.0, 34.0, 38.0, 42.0, 46.0, 50.0, 54.0, 58.0};
    Float_t dtimes[12] = {0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5};

     // Create timegraph
    TGraphErrors * gr = new TGraphErrors(12,peaks,times,dpeaks,dtimes);
    gr->SetMarkerColor(kRed);
    gr->SetMarkerStyle(20);
    gr->SetMarkerSize(0.7);

    TF1 tFunc("tFunc", lin, 500, 3000, 2);
    tFunc.SetParNames("Slope", "Offset");
    tFunc.SetParameters(0.02,0);
    tFunc.SetLineColor(kBlue);
    tFunc.SetLineWidth(2);
   
    gr->Fit(&tFunc);



    // Setze Styles
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(0);

    // Setze Beschriftungen und Styles
    gr->SetTitle("Time Calibration");
    gr->GetXaxis()->SetTitle("channel");
    gr->GetYaxis()->SetTitle("time / ns");

    // Erzeuge das Fenster
    TCanvas* Canva = gROOT->Find Object("Canva");
    if(!Canva) Canva = new TCanvas("Canva", "Versuch_05 TimeCal", 1024, 768);
    else Canva->Clear();
    Canva->cd(1); // Waehle das Feld an


    // Male Graph
    gr->Draw("APE");
    tFunc.Draw("Same");
}





/// Define fit function for Cobalt-60 angle correlation
double co60ang(double* vars, double* pars) {
  return pars[1]*TMath::Cos(vars[0])*TMath::Cos(vars[0]) + pars[2]*TMath::Cos(vars[0])*TMath::Cos(vars[0])*TMath::Cos(vars[0])*TMath::Cos(vars[0]) + pars[0];
}


/// Read angle correlation data from file
void angleCor(char* filename, int const N)
{

  // Define constants
  double const Tc = 10e-9;   /* [s] Coincidence Time       */
  double const dTc = 2e-9;   /* [s] Coincidence Time Error */
  double const T = 1800.0;   /* [s] Measurement Time       */
  double const dT = 1.0;     /* [s] Measurement Time Error */
  double const dang = 0.5;   /* [°] Angle Error            */

  // Declare local fields
  Double_t angle[N], N1[N], N2[N], Nk[N];
  Double_t dangle[N], dN1[N], dN2[N], dNk[N];

  Double_t Z1[N], Z2[N], Zk[N], Ze[N], Zz[N];
  Double_t dZ1[N], dZ2[N], dZk[N], dZe[N], dZz[N];

  FILE * outfile = fopen("Co60_ang.tex","w");
  FILE * outfile2 = fopen("Co60_ang.val", "w");
  FILE * outfile3 = fopen("Co60_ang.nrm", "w");

  // File Opening Reading and Closing
  FILE* file = fopen(filename, "r");
  for(Int_t i=0; i < N; i++){
    fscanf(file, "%lf %lf %lf %lf", &angle[i],&N1[i],&N2[i],&Nk[i]);
  }
  fclose(file);

  // Data processing
  for(Int_t i=0; i < N; i++){
    angle[i] = angle[i]*TMath::DegToRad();

    // Calculate erros
    dangle[i] = dang*TMath::DegToRad();
    dN1[i] = N1[i]/TMath::Sqrt(N1[i]);
    dN2[i] = N2[i]/TMath::Sqrt(N2[i]);
    dNk[i] = Nk[i]/TMath::Sqrt(Nk[i]);

    // Calculater count rates;
    Z1[i] = N1[i] / T;
    Z2[i] = N2[i] / T;
    Zk[i] = Nk[i] / T;
    dZ1[i] = TMath::Sqrt( dN1[i]*dN1[i]/(T*T) + N1[i]*N1[i]*dT*dT/(T*T*T*T) );
    dZ2[i] = TMath::Sqrt( dN2[i]*dN2[i]/(T*T) + N2[i]*N2[i]*dT*dT/(T*T*T*T) );
    dZk[i] = TMath::Sqrt( dNk[i]*dNk[i]/(T*T) + Nk[i]*Nk[i]*dT*dT/(T*T*T*T) );
    
    // Calculate true coincidences
    Ze[i] = Zk[i] - Tc*Z1[i]*Z2[i];
    dZe[i] = TMath::Sqrt( dZk[i]*dZk[i] + Z1[i]*Z1[i]*Z2[i]*Z2[i]*dTc*dTc + Tc*Tc*Z2[i]*Z2[i]*dZ1[i]*dZ1[i] + Tc*Tc*Z1[i]*Z1[i]*dZ2[i]*dZ2[i] );

  }
  
  // Norm to Ze(90°)
  // Cobalt-60 only
  for(Int_t i=0; i < N; i++) {
    dZe[i] = TMath::Sqrt( dZe[i]*dZe[i]/(Ze[4]*Ze[4]) + Ze[i]*Ze[i]*dZe[4]*dZe[4]/(Ze[4]*Ze[4]*Ze[4]*Ze[4]) );
    Ze[i] = Ze[i]/Ze[4];
  }
  //dZe[4] = dZe[4]/Ze[4];
  //Ze[4] = 1.0;

  // Write Latex File
  for(Int_t i=0; i < N; i++) {
    fprintf (outfile, "%.1f & %d & %d & %d \\\\ \n",angle[i]*TMath::RadToDeg(),N1[i],N2[i],Nk[i]); 
    fprintf (outfile2, "%.1f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\ \n",angle[i]*TMath::RadToDeg(),Z1[i],dZ1[i],Z2[i],dZ2[i],Zk[i],dZk[i]);
    fprintf (outfile3, "%.1f & %.2f & %.2f \\\\ \n", angle[i], Ze[i], dZe[i]);
  }
  fclose(outfile);
  fclose(outfile2);
  fclose(outfile3);

  // Create TGraphErrors representation
  TGraphErrors* graph1 = new TGraphErrors(N,angle,Ze,dangle,dZe);
  graph1->SetTitle("True Coincidences");
  graph1->GetXaxis()->SetTitle("#Theta / rad");
  graph1->GetYaxis()->SetTitle("W(#Theta)");
  graph1->SetMarkerStyle(2);
  graph1->SetMarkerSize(1);

  // Create Fit-Function for Natrium-22
  TF1 fit1("fit1",gaussian,145,195,4);
  fit1.SetParNames("Norm", "Mean", "StdDev", "Offset");
  fit1.SetParameters(35,180,10,0);
  fit1.SetLineColor(kRed);

  // Create Fit-Function for Cobalt-60
  TF1 fit2("fit2", co60ang, 1.5, 3.2, 3);
  fit2.SetParNames("Offset", "a", "b");
  fit2.SetParameters(0.97,0.125,0.4167);
  fit2.SetLineStyle(1);
  fit2.SetLineColor(kRed);

  // Create theorectial function for Cobalt-60
  TF1 theo("theo", "1.0+cos(x)*cos(x)/8+cos(x)*cos(x)*cos(x)*cos(x)/24",1.5,3.2);
  theo.SetLineColor(kBlue);
  theo.SetLineStyle(2);
  theo.SetLineWidth(1);

  // Perform Fit
  //graph1->Fit(&fit1,"R");     /* Natrium-22 */
  graph1->Fit(&fit2,"R");     /* Cobalt-60 */

  // Setze Styles
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  // Create the canvas for painting
  TCanvas* c1 = gROOT->Find Object("c1");
  if(!c1) c1 = new TCanvas("c1", "Versuch_05", 1024, 768);
  else c1->Clear();
  c1->cd(1);

  // Paint the graph
  graph1->Draw("APE");
  theo.DrawClone("Same");
}  


/// Read data from file and plot a histogram
void histo2(char* filename)
{
    int const N = 4096;  

    // File Opening Reading and Closing
    FILE* file = fopen(filename, "r");
    char dummy[12];
    // Lese Daten
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
	if (channel[i] > 2550) value[i] = 0;
    }
    fclose(file);

    // Create 1-D Histogram
    TH1F* Hist = gROOT->FindObject("Hist");
    if(!Hist) Hist = new TH1F("Hist", "Histogramm", N, 0, N);
    else Hist->Reset("M");

    // Fill histogramm
    for(Int_t i=0; i < N; i++){
      Hist->Fill(channel[i],value[i]);
    }

    Hist->Rebin(2);

    // Setze Styles
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(0);

    // Setze Beschriftungen und Styles
    Hist->SetTitle("ADC Spectrum");
    Hist->GetXaxis()->SetTitle("channel");
    Hist->GetYaxis()->SetTitle("counts");

    // Erzeuge das Fenster
    TCanvas* Canva = gROOT->Find Object("Canva");
    if(!Canva) Canva = new TCanvas("Canva", "Versuch_05", 1024, 768);
    else Canva->Clear();
    Canva->cd(1); // Waehle das Feld an

    // Male Graph
    Hist->Draw("PL");
}
